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Ab-initio simulations of quantum transport commonly focus on a central region which is consid- 
ered to be connected to infinite, periodic leads through which the current flows. The electronic 
structure of these distant leads is normally obtained from an equilibrium calculation, ignoring the 
self-consistent response of the leads to the current. We examine the consequences of this, and 
show that the electrostatic potential, A<j>, is effectively being approximated by the difference be- 
tween electrochemical potentials, Afi, and that this approximation is incompatible with asymptotic 
charge neutrality. In a test calculation for a simple metal-vacuum-metal junction, we find large 
errors in the non-equilibrium properties calculated with this approximation, in the limit of small 
vacuum gaps. We provide a simple scheme by which these errors may be corrected. 

PACS numbers: 05.60.Gg, 73.40.-c, 41.20 



I. INTRODUCTION 

The last fifteen years have seen considerable progress 
in the simulation of non-equilibrium many-electron 
nanoscale systemsiaSiMi&i&i typically using the Lo- 
cal Density Approximation (LDA) to Density Func- 
tional Theory (DFT)^ within the framework of non- 
equilibrium Green's functions (NEGF)^. In these simu- 
lations the systems under study consist of two electrodes 
(or leads), placed to the left and to the right of an ac- 
tive central region which contains a molecule and parts 
of the left and right electrode. Starting from the un- 
connected electrode-central region-electrode system with 
each of the electrodes itself in equilibrium but not in 
equilibrium with each other— the NEGF formalism then 
provides the formal apparatus to switch on the contacting 
terms of the Hamiltonian adiabatically, causing a current 
to flow through the system. Associated with this current 
and related to the resistance of the molecule there is a 
"resistivity dipole" arising from the newly induced charge 
density, that causes the electrostatic potential to drop in 
the neighbourhood of the central region. The magnitude 
of the drop in the self-consistent electrostatic potential 
is essentially fixed by a charge neutrality condition, i.e., 
the fact that the asymptotic electrode regions must them- 
selves be charge-neutral since a net charge would cause 
the electrostatic potential to diverge. In the case of jel- 
lium electrodes, this charge neutrality condition acquires 
a strict local form 2 ^ since, asymptotically, the electron 
density exactly cancels the background density at any 
point. 

It should be pointed out that most practical implemen- 
tations of the NEGF formalism cannot, strictly speak- 
ing, properly take into account the drop in the electro- 
static potential between the leads, since the adiabatic 
switching of the contacting terms of the Hamiltonian (i.e. 
the perturbation) changes the density and the response 
to that change cannot be described inside the realm of 
static DFT. Furthermore the lead self-energies, which 



describe the coupling between the leads and the cen- 
tral region, are commonly obtained from an equilibrium 
calculation^. At the level of the Hartree approximation 
and in the non-equilibrium regime, this means that the 
leads do not respond to the flow of charge induced by 
the applied bias voltage; the electrochemical potentials 
remain fixed to their equilibrium values and the drop 
in the self-consistent electrostatic potential, A0, is, as 
discussed below, effectively being approximated by the 
difference between electrochemical potentials, A/ji*^*i£. 
At the level of the Hartree-Fock approximation this lack 
of asymptotic self-consistency also implies that the non- 
equilibrium Fock operator deep inside the leads would be 
equal to the equilibrium one, which is clearly not the case 
since the Fock operator depends on the non-equilibrium 
occupancies of the current carrying states, which are dif- 
ferent from the equilibrium ones everywhere. However, 
we will not further discuss the effects of the lack of asymp- 
totic self-consistency in the Fock operator. The most se- 
vere effects appear already in the self-consistent Hartree 
potential which we will use as an illustration in our pa- 
per. We would like to note that non-partitioned NEGF 
approaches, as suggested by Cinii£ and late elaborated 
by Stefanucci and Almbladhi^i^, are in theory free from 
these objections as they focus on the evaluation of the 
non-equilibrium Green's function in the whole transport- 
ing system. 

The relation between the electrostatic drop and asymp- 
totic charge neutrality is already implicit in the original 
form of the Landauer formulaiSiii / = ^A0, and has 
been explored by some authors over the yearsi&iS until 
very recentl}^. In Ref. |2fJ the authors further clarify 
the distinction between the difference between the elec- 
trochemical potentials of the left and right electrodes, 
Afi, and the drop in the electrostatic potential, Acj), as 
well as the role played by the geometry. However they 
do not discuss in detail the validity of the approximation 
A^ = Ac/). 

To fix ideas consider the biased system depicted in 
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FIG. 1: (color online). Schematic illustration of the self-consistent electrostatic potential (thick solid line) together with its drop 
A<f>, and the difference between electrochemical potentials A/i = pa — pl- In (a) and (b) Ap = Acf>. Solid (dashed) arrows are 
indicative of the transmission and reflection amplitudes of right-going (left-going) states. In (a) a thick barrier with neglible 
transmission between pl and pa is shown. In this case for the system to be neutral in its far left and right regions Ap = A(j>. 
The thinner barrier in (b) allows for larger values of the transmission coefficient and hence, since there is a net and significant 
flow of particles from left to right, there local density of states is smaller on the left than on the right. If A/i = A<f> this case 
leads to a net accumulation of charge on the right region and a depletion of it in the left region. Therefore A<f> has to deviate 
from A/i in the way schematically shown in (c). 



Fig. ^ which is translationally invariant in two of the 
three spatial directions and has some localized inhomo- 
geneity along the z-direction. Since left- and right-going 
scattering states are occupied up to two different electro- 
chemical potentials, fiL and p,R respectively, a current 
flows along z, and therefore there is an electrostatic po- 
tential drop A<p. If we assume that A/i = (Xr—[1l = A<f>, 
then it is easy to show that in the asymptotic electrode 
regions the electronic density, p, is in the case of jellium 
electrodes, given bjii2ii2i2i: 

p(z — > ±oo) = pb(z — > ±oo) 

/•+0O 

± / dk ± T(E z )[f R (E z ) - f L (E z )](l) 
Jo 

where pb{z) is the background density, k± are the mag- 
nitudes of the z-component of the momentum in the 
asymptotic right (+) and left (— ) electrode regions, re- 
spectively; E z (k±) = k 2 _/2 = fc 2 /2 — A<j) is the energy 
associated with the motion in the direction of the cur- 
rent; T[E Z ) is the usual transmission probabilities and 
fa and /l are the equilibrium Fermi-Dirac distributions 
for right- and left-going electrons averaged over the com- 
ponents of the momentum perpendicular to the direction 
of the current, each of which is characterized by an elec- 
trochemical potential /j-r(l)- The latter, as functions of 
fe_ , are given bjii&iiLSI: 

f L , R (k-) = J^(*£,r - *£)e(fc_ - k L<R ), (2) 

where fc_L._R are the asymptotic Fermi wave-vectors for 
left- and right-going electrons respectively, in the left 
asymptotic region. Therefore, when imposing A/i = A<p, 
associated with the presence of a current from, say, left 
to right, there is a charge depletion in the asymptotic 
left electrode region and a charge accumulation in the 
asymptotic region of the right electrodsi2ii2i22i2i. There- 
fore it is clear that the drop in the self-consistent electro- 
static potential is necessarily different from A/i. It is then 



surprising that in many state-of-the-art ab-initio quan- 
tum transport simulations the approximation A/i = Acj) 
is used without further explanation or comments on its 
validityi^ii^. Of the few that have considered this prob- 
lem let us mention P6tai&, who introduces a drift in the 
electronic distribution functions of left- and right-going 
electrons so that the asymptotic electrode regions remain 
neutral, Bokes and Godbyi^ whose proposed method is 
applied in this work and Di Ventra and Lang 22 who renor- 
malize the electron densities deep in the jellium elec- 
trodes. 

From the above given discussion it should be clear that 
asymptotic self-consistency (i.e. the full non-equilibrium 
Green's function of the leads) is generally needed in or- 
der to describe the drop in the electrostatic potential 
correctly. The lack of asymptotic self-consistency natu- 
rally leads to the approximation A/i = Aip which is in- 
compatible with charge neutrality in the asymptotic lead 
regions. The rest of this paper is organized as follows: 
In Section II we show explicitly the effect that the lack 
of asymptotic self-consistency has in the calculated non- 
equilibrium properties of a simple jellium metal- vacuum- 
metal junction. In Section III we express these ideas in 
the lenguage of NEGF's, showing which terms are com- 
monly neglected and providing a simple recipe to incorpo- 
rate them for the case of jellium electrodes. We conclude 
in Section IV. 



II. EXAMPLE: ASYMPTOTIC CHARGE 
NEUTRALITY IN A JELLIUM 
METAL- VACUUM-METAL INTERFACE 

Neglecting the non-equilibrium contributions to the 
Green's function of the leads results, at the level of the 
Hartree or DFT-LDA approximations, in the approxi- 
mation A/i = A(f>, which in turn is incompatible with 
asymptotic charge neutrality in the leads. We now show 
the effects of this approximation in the calculated prop- 
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erties of a non-equilibrium metal- vacuum-metal junction. 
Even though we will use a scattering-state-based ap- 
proach its equivalence with the NEGF formalism may 
be noted. 

The jellium model of the mctal-vacuum-metal 
interfacei*2iSi is defined in terms of the background den- 
sity: 

PB (z) = no[9(-z) + 9(z - L)}, (3) 

where n = and L is the length of the vacuum 

gap. For this system we solve the Kohn-Sham equations 
self-consistently using the Perdew-Zunger— parametriza- 
tion of the LDA exchange-correlation potential. His- 
torically this system was the first to be studied using 
conventional ab-initio techniques in a non-equilibrium 
regime 1-24 , and, for our purposes, constitutes a simple 
system for which the electrostatic effects under study 
arise in the most transparent manner. 

In order to ensure charge neutrality in the asymptotic 
electrode regions, for a given value of A<p we need to find 
kn and fc^ in Eqs. @ such that 

n = p(z —> ±00), (4) 

are satisfied, kn and fej, are related to the electrode elec- 
trochemical potentials simply by 



Therefore, at each step of the self-consistency cycle, we 
solve the Poisson equation with Dirichlct boundary con- 
ditions, fixing A(j) and calculating the corresponding A/i 
that ensures that the asymptotic left and right electrode 
regions remain neutral. Strictly speaking this procedure 
is only justified in the case (such as the metal-vacuum- 
metal junction) that there is a one-to-one correspondence 
between A/i and A(j>, i.e., there is a one-to-one corre- 
spondence between the applied bias and the current. For 
this particular case our method is equivalent to the al- 
ternative one of fixing A/i and calculating Acj) by solv- 
ing the Poisson equation with von-Neumann boundary 
conditions^. 

When studying the influence of asymptotic charge 
neutrality in the calculated non-equilibrium properties 
we solve the Kohn-Sham equations self-consistently us- 
ing the procedure described by McCann and Brown 1 , 
with the approximation A/i = Acf> and without it, us- 
ing the asymptotic neutrality condition Eq. lfl|. In par- 
ticular the Poisson equation is solved using Nieminen's 
methodiiS^, which greatly stabilizes the iterative pro- 
cess and fastens the convergence. 

We next present a set of results for a symmetric metal- 
vacuum-metal junction with electrodes characterized by 
r s = 4 2 », focusing in the differences between A</> and Ap 
as a function of the electrode-electrode distance and the 
effect that the lack asymptotic neutrality has in the cal- 
culated resistivity dipoles and current densities. Fig. |3 
shows a linear relation between Ap and A<f> for different 




A<t>/E, 



FIG. 2: (color online) The figure shows the difference between 
electrochemical potentials Ap as a function of the drop in the 
electrostatic potential A(j>, in units of the equilibrium Fermi 
energy, for different values of the electrode-electrode distance, 
L. L — 0.25r s (solid line); L — 0.5r s (dashed); L — lr s (short 
dashes); L = 1.5r s (dots); L — 2r s (dot-dashed). A/i ~ 
A(f> only for large electrode-electrode spacings. For reference 
E F = 3.131 eV 

lengths of the vacuum gap. All the lines fall between 
two limiting ones: A<p = for L = and Aip — Ap for 
L — > 00 as expected. For L = the system is homoge- 
neous along the z direction and hence there is no electro- 
static drop. As the distance between electrodes increases 
the transmission coefficient T(E Z ) decreases, the current 
decreases and p(z — > ±00) — > no, therefore the electro- 
static drop and the difference between electrochemical 
potentials are approximately equal in that limit. In fact 
A/i w A(f) for L > 2r s . However one should note that 
there are no molecular conducting channels present in 
our model. If these were present and open, then the de- 
viation from A/i = A<j) should, according to Eq. QJ, be 
larger at a fixed value of the electrode-electrode separa- 
tion. Fig. shows calculated resistivity dipoles defined 
as: 

6p(z)=p(z,A<j>^0)-p(z,A<p = 0), (6) 

for different values of the applied bias. In Fig |3Ja) 
the dipoles were calculated within the approximation 
A/i = Acf>. Enforcing this boundary condition when solv- 
ing the Poisson equation leads to the appearance of un- 
physical charges which are placed at the edges of the nu- 
merical grid used in the calculations. These spurious con- 
tributions to the induced density disappear as A/i — > 
or L — > 00, since in these limits A/i = A<f>. In Fig 0{b) 
we show the calculated resistivity dipoles by choosing 
and kn so that Eqs. 0] are satisfied. The induced density 
goes smoothly to zero as z — > ±00 even at small values of 
L and relatively large values of A</>. Finally in Fig. 0] we 
present calculated J — AV curves (with A V equal to A/i 
or A(j) depending on the pair of curves being compared) 
for different electrode-electrode spacings. Large differ- 
ence between the J — A(f> — Ap curve and any of the 
other two are present at small electrode-electrode sepa- 
rations. As argued above, as the separation between the 
electrodes becomes larger all three curves converge into 
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FIG. 3: (color online) (a) resistivity dipoles calculated for dif- 
ferent values of A(p with L = 1.5r s , using A<f> ~ Afi. The ar- 
rows indicate the presence of unphysical charges at the edges 
of the numerical grid; (b) same as in (a) but calculated using 
our neutrality scheme. Solid line A<P/Ef = 0.075, dashed 
A<f>/E F = 0.05, dotted A(j>/E F = 0.025. The vertical lines 
indicate the positions of the edges of both jellium surfaces 

a single one. Table [I] contains numerical values of the ra- 
tios between linear response conductances (calculated at 
a small but finite bias) with and without the approxima- 
tion A/i = A<p. The performance of this approximation is 
poorer for larger values of the conductance, as expected. 



TABLE I: Ratios between the calculated linear response con- 
ductances per unit area Gip = J/Afi, Gap = J/A<fi and 
GAfi=A<p = J/ AV with AV — Afj, = A<p, together with their 
corresponding value of the transmission coefficient evaluated 
at E z = E F . 



L/r s 
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1.0 


0.89 
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1.00 


1.00 
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III. RELATION TO THE PARTITIONED NEGF 
APPROACH 

Even though the problem of asymptotic charge neu- 
trality for jellium electrodes is clear and directly solv- 
able through the scattering states as it is done in this 





paper, the issue becomes more involved within the par- 
titioned NEGF approach. We proceed to show explic- 
itly the terms that are usually neglected in practical im- 
plementations and to provide a simple scheme by which 
asymptotic self-consistency may be easily implemented. 

Let us consider a partitioned system's Hamiltonian 
into two semi-infinite leads ( left (L) and right (R) ) and 
finite central region (C). The Green's functions we con- 
sider below are all defined on the Keldysh contour with 
the complex time- variable t— and are meant to represent 
matrices with two indexes m, n belonging to some com- 
plete set of spatially localised basis functions. Each of 
these can belong to any of the above introduced regions 
L,R or C. We employ the notation 

[G a ]„ m = -z(T T {c„ eQ (r)cJ riea (r')}) a = L,R,C (7) 

for the disconnected systems and 

[G a0 ] nm = ~i(T T {c n<Ea (r)cl ef3 (r')}) a,/3 = L, R, (7(8) 

for the contacted ones. We also use 1 = S nm 6(T — r'). 

Before we turn on the contacting between L and C 
and C and R the uncontacted Green's functions fulfill 
the equations of motion 

(id T -H L )G L = 1 (9) 
(id T -H c )G c = 1 (10) 
(id r -H R )G R = 1. (11) 

This can be also written in a block form as 



f H L \ 
id T l - H c 

\ H R J 

or more concisely as 



G L \ 
G c = 1(12) 
G R J 



{id T - H)G° = 1. (13) 

Next we turn on the interaction terms which couple the 
left and central and the right and central parts, written 
as Vl and V R respectively. The coupling, however, in- 
duces also a change in Hamiltonians via the change of the 
density in the Hartree and the exchange-correlation po- 
tentials which we together write as 5H a for a = L,C,R. 
Using the block Green's functions we have 

(id T - H - SH - V)G = 1. (14) 

where 

/H° L +6H L V L \ 

H + SH = Vl H c + SH C V R 

V Vl H R + SH R J 

and 

(Gll Glc Glr \ 
Gcl Gcc Gcr J ■ 
Grl Grc Grr I 
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FIG. 4: (color online) J — AV characteristics for different electrode-electrode spacings. (a) L = lr a ; (b) L = 1.5r s ; (c) L = 3r s . 
For the solid lines AV = A</> = A/x, for the dotted lines AV = A0 and for the dashed lines AV = A/x 



The solution of the Eq. (|14f> can be written using Eq. 
in the form of the Dyson equation as 

G = G° + G° [SH + V]-G, (15) 

where • stands for the integral along the Keldysh contour 
over an internal time variable. 

The Green's function of the central region corresponds 
to the finite system and is usually solved numerically in 
the self-consistent manner. It solves the Dyson equation 

Gcc = Gc + GcSHc ■ Gcc 

+ GcVl-Glc + GcVr-Grc (16) 

which can be found as the 'CC component of Eq. (|15|l . 
To have a closed system of equations we need to find 
Glc and Grc- These are similarly given as 'LC and 
'RG components of Eq. H15[l as 

Glc = GlVl ■ Gcc + GlSHl ■ Glc, (17) 

and 

Grc = GrVr ■ Gcc + GrSHr ■ Grc- (18) 

The latter two can be formally inverted (in m, n as well 
as in r) to give 

G LC = [1-Gl6HlY 1 -GlVl-Gcc (19) 
Grc = [1-GrSHr]- 1 -GrVr-Gcc- (20) 

Combining Eqs. ((Tfj|l . (|T§|) and lj2U|) we finally obtain 

Gcc — Gc + Gc {SHc + ~£>cc) ■ Gcc, (21) 
Sec = V L [1 - GlSHl}- 1 ■ G L V L 

+ VrII-GrSHr]- 1 -GrVr (22) 

= VlGlVl + VrGrVr, (23) 

where we have introduced the self-energy £cc represent- 
ing the leads for the central region and defined new aux- 
iliary Green's functions Gl /r that fulfill the equations of 
motion 

[id T - H a - 5H a ) G a = l, a = L, R. (24) 



The G a need to be such that they correspond to the equi- 
librium situation with 5Hl/r being turn on but with 
chemical potentials (or Fermi-Dirac occupation factors 
/l and /r) being kept the same as in Gl/r- For this 
reason the only significance of the auxiliary Green's func- 
tions is their presence in the expression for the self-energy 
'Sec, having no other direct physical meaningSi. Using 
the calculated Gcc in the Eq. I|19|) one can now employ 
the usual derivation of the expression for the current in 
terms of G^^^i 

I=^{Tt[VG< c ]}. (25) 

The usual treatments'^ do not consider the change 
in the leads' Hamiltonians due to the change in density 
SHl/r that arises in the non-equilibrium regime. This 
results in a simplification of our equations since the re- 
solvent operators [1 — GrSHr]' 1 and [1 — GlSHl] 1 do 
not have to be calculated in the non-equilibrium regime. 
We note that the evaluation of these would require a 
complete calculation of Gll and Grr for the contacted 
system since these give the non-equilibrium density that 
in turn determines the changes SHl/r- This can be 
most easily achieved for semi-infinite jellium electrodes 
which is equivalent to the calculations presented in the 
previous section of this paper (see also reference^ for 
model ID cases treated directly using the non-partitioned 
NEGF formalism). However completely ignoring these 
changes results in a violation of asymptotic charge neu- 
trality in the leads! A practically useful scheme would 
be to consider this change to be just a constant shift, i.e. 
5H l /r = SUl/r- This is exact for jellium electrodes if 
leads are taken to be sufficiently far from the constriction. 
As it can be seen from Eq. (j24(l , this shift moves the bot- 
toms of the bands in the leads' density of states, which 
is the only characteristic of leads that eventually enters 
into the final expression for the current Eq. (|25|l . The 
relation between the drop in the chemical potentials, the 
drop in the induced potential and these shifts is simply 

A0 + SU L - SUr = fi R -iil, (26) 

which is also clear from Fig. lc. Unfortunately this rela- 
tion fixes only the differences between the shifts SUl/r- 
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We still need to have one more equation to determine 
them uniquely. This would come from the computation 
of non-equilibrium local charge neutrality in one of the 
electrodes. 



IV. CONCLUSIONS 

We have shown that in an exact time-dependent 
density-functional formulation of the partitioned 
Keldysh-NEGF approach the changes in the Hamil- 
tonian of the leads due to the contacting needs to be 
included. It is important for a correct description of the 
electrostatic potential profile at large currents or junction 
with transmission close to one. Using a simple jcllium 
model of a biased metal-vacuum-metal junction we have 
examined quantitatively the effects of fixing A/z = A(f>. 
Significant differences between the non-equilibrium 



properties calculated using this approximation and a 
more reasonable treatment of the electrostatic problem 
based on asymptotic charge neutrality arise in the 
limit of small electrode-electrode separation. These 
effects would be even more pronounced for resonant 
molecular junctions where both electrostatics and high 
conductance acting simultaneously may significantly 
influence the I — V characteristics of the system. 
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